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We carry out a numerical study of the bi-partite entanglement entropy in the gapped regime of 
two paradigmatic quantum spin chain models: the Ising chain in an external magnetic field and 
the anti-ferromagnetic XXZ model. The universal scaling limit of these models is described by 
the massive Ising field theory and the S'C/(2)-Thirring (sine-Gordon) model, respectively. We may 
therefore exploit quantum field theoretical results to predict the behaviour of the entropy. We 
numerically confirm that, in the scaling limit, corrections to the saturation of the entropy at large 
region size are proportional to Ko{2mr) where m is a mass scale (the inverse correlation length) 
and r the length of the region under consideration. The proportionality constant is simply related 
to the number of particle types in the universal spectrum. This was originally predicted in [1, 2] 
for two-dimensional quantum field theories. Away from the universal region our numerics suggest 
an entropic behaviour following quite closely the quantum field theory prediction, except for extra 
dependencies on the correlation length. 
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Introduction and discussion. Entanglement is a 
fundamental property of the state of a quantum system. 
In the context of quantum computation, it is a crucial 
resource [3]. Conceptually, it characterizes the structure 
of quantum fluctuations in a more universal way than 
other widely studied objects such as correlation func- 
tions. Developing theoretical measures of entanglement 
is therefore important to further understand the struc- 
ture of quantum states, a problem of particular inter- 
est and difficulty for quantum many-body systems. A 
popular measure of entanglement is the bi-partite entan- 
glement entropy. It measures the entanglement between 
two complementary sets of observables in a quantum sys- 
tem [4] . Interestingly, the entanglement entropy exhibits 
universal behaviour near quantum critical points: it has 
features which do not depend on the details of the model 
but rather on its universality class. In the last decade, 
this property has made the study of the entanglement 
entropy a very active field of research. 

A fertile testing ground for these ideas is the study of 
quantum spin chains. A quantum spin chain is a one- 
dimensional array of particles with spin degrees of free- 
dom and, for our purposes, with nearest- (or few-nearest- 
) neighbour interactions (this is a local spin chain). 
Quantum spin chains have been realized in experiments 
[5]. They provide ideal toy models for the study of a 
whole range of physical phenomena, as they describe in- 
teracting many-body systems yet are simple enough to 
allow for the computation of many quantities. This is 
especially true for integrable spin chains, since integra- 
bility gives rise to the complete characterization of the 
energy spectrum and states through techniques such as 
the coordinate and algebraic Bethe ansatz (see e.g. [6]). 
Quantum spin chains are extensively studied in the con- 
text of quantum information science [7] and their bipar- 
tite entanglement for blocks of consecutive spins has been 
investigated in many works (see e.g. [8-14]) both numer- 
ically and analytically. Most of these works concentrate 
on exact critical points. 



At the same time, in 1+1-dimensional quantum field 
theory (QFT) several universal results have been ob- 
tained for the entanglement entropy. Two of these results 
deserve special attention. First, the logarithmic growth 
of entanglement as the block size increases in conformal 
field theory (CFT), which is controlled by the central 
charge [15, 16] . This behaviour has been verified for many 
critical quantum chains [8-11, 14]. Second, the approach 
to saturation in massive QFT, which is controlled by the 
mass spectrum [1, 2]. This has not been observed yet in 
quantum spin chains. 

In the present work we investigate the entanglement 
entropy in two paradigmatic spin chain models (the Ising 
model and the XXZ model) near to, but not at, critical- 
ity. We numerically analyze for the first time the uni- 
versal behaviour of the entanglement entropy of a block 
of length L in the near-critical scaling limit. This scal- 
ing limit is expected to be described by massive QFT. In 
particular, we confirm massive QFT predictions [1, 2]. 

The entanglement entropy is the von Neumann entropy 
of the reduced density matrix of a state |^) with respect 
to a tensor factor of the Hilbert space TL ^ A<S^ B: 



S 



-Tta{pa log pa) with pA-TrB|*)(*|. (1) 



Let 1^) be the ground state of a quantum spin chain. 
In general, it is expected that, thanks to locality, the 
entanglement entropy of a continuous block A of length 
L saturates as L — >■ oo: it only receives contributions 
from entanglement between spins surrounding boundary 
points of A. But this is not so at second order phase tran- 
sitions. Assume that the Hamiltonian is parametrized by 
h, and that there exists a value he corresponding to a 
quantum second order phase transition (critical point). 
These critical points are particularly interesting, as they 
characterise collective quantum behaviours. At h — he, 
and in the thermodynamic limit where the length of the 
chain A^ — ;► oo, the correlation length ^ of the ground 
state is infinite, and the gap oc ^~^ between the ground 
state and the continuum of excited states vanishes: the 



system is critical. If the dynamical exponent is z = 1 (so 
that the dispersion relation is linear), then the macro- 
scopic, low-energy properties are universal and described 
by a CFT, see FIG. 1. 

For a block A of length L, the entanglement entropy 
S diverges for large L as [15, 16]: 
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where c is the central charge of the CFT, e is a short- 
distance scale (here and below taken as the inter-site 
spacing) and C2 is a non-universal constant. The diver- 
gence occurs because with ^ = oo sites far apart are 
entangled. Thus we may extract the value of the central 
charge by studying the divergence of the entanglement 
entropy in critical spin chains. 

In obtaining (2), 
the order of limits, 
first h —i' he then 
L ^ e, is very im- 
portant. Instead, let 
us take h ^ he and 
L ^ e "simultane- 
ously": with L/^ =: 
mr fixed. This is 
the near-critical scal- 
ing limit, the shaded 
region depicted in 
FIG. 1. It is uni- 
versal, described by 
a massive QFT corre- 
sponding to a pertur- 
bation off the orig- 
inal critical point. 
The mass scale m 
may be identified with ^~^. In this limit, the gap of 
the spin chain is infinitesimally small as compared to mi- 
croscopic energies (like the inter-site interaction energy), 
but the observation length is so large as to probe only the 
low-energy, universal, collective excitations of the quan- 
tum chain just above the gap; hence the gap still has 
an effect. The entanglement entropy then has the form 
of a "saturation term" (c/3) log(^/e) -I- ci which diverges 
in the scaling limit, where ci is another non-universal 
constant, plus a universal scaling function f(mr). The 
approach to saturation as mr — > oo is given by an expo- 
nential decay of f{mr) which is solely determined by the 
spectrum of masses {m;} of the QFT [1, 2]: 



FIG. 1: Quantum field theory re- 
produces the universal scaling limit 
of physical quantities: observation 
energies and the gap between the 
ground state and the first excited 
state are both very small as com- 
pared to the microscopic energy 
scales, like inter-site interaction en- 
ergies. 
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-ci-l^A-o(2rmO + 0(e-3™) (3) 



where m = mi is the mass of the lightest particle. Here 
Ko{z) is the modified Bessel function. Thus, in the near- 
critical scaling limit, the entanglement entropy of the 
chain encapsulates information about the mass spectrum 
of the QFT. Further, despite the non-universality of ci 
and C2, their difference is universal: 



In massive QFT, U is related to the expectation value of 
a branch point twist field, as described in [1]. In contrast 
to the large number of works dealing with critical spin 
chains, the absence of studies of the near-critical scaling 
regime precludes a comparison with massive QFT results. 
The present work intends to fill this gap. 



Ising model. This model has Hamiltonian 



T ^ 



haf 



(5) 



where af are Pauli matrices acting on site i and we con- 
sider periodic boundary conditions <7i+N = <^i- The "mi- 
croscopic" energy scale indicated in FIG. 1 is the coupling 
constant J > 0, and h plays the role of an external mag- 
netic field. It is well-know (see e.g. [17]) that this model 
has a quantum critical point for h = he = I described 
by a CFT with c = ^. The exact correlation length of 
the chain is given by ^^^ = log ft- [18] and, as expected, 
diverges at the critical point. For h > 1 the system is 
gapped and the ground state is unique (no symmetry 
is broken). Taking the scaling limit described above, a 
massive QFT is obtained: the free massive relativistic 
Majorana fermion. 

A numerical study of the Ising model has the advan- 
tage that a free fermion map can be used to perform 
computations in the thermodynamic regime. This tech- 
nique is used and explained in much detail in [9] and it 
is based on the use of Toeplitz determinants [11]. We 
can then carry out extremely precise numerics on an ex- 
actly infinite chain for very large values of L and very 
small values of £,~^, which are numerically inaccessible 
for more complex models. 

A numerical computation at the conformal critical 
point yields the expected behaviour (2). For the Ising 
model, the exact value oi U = —0.131984... (we omit 
here the exact expression) was first evaluated in [1]. With 
ci = log 2/2 calculated in [19] this gave, using (4), the 
constant C2 = 0.478558.... This was found to agree ex- 
tremely well with the numerical results obtained in [9]. 
Our numerics also confirm this value, see FIG. 2. 
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FIG. 2: (color online) Entanglement entropy of the Ising 
model for several values of h. At the critical point h = 1, 
the entropy is very well described by the function (2) with 
c = 0.500003 and C2 = 0.478551 (red dots). Away from the 
critical point we observe the rapid saturation of the entropy 
to a constant value. 



Away from criticality, the spectrum of the QFT has a 
single particle. Hence formula (3), with only one term 
in the sum, should describe the large-Tnr behavior in the 
limit h — >■ 1+. In our numerics, we take finite, increasing 
values of ^, and each time the large-L behaviour is fitted 
to 

5=iloge + ci(0-^A'o(2L/0 (6) 

(we choose e = 1 and use mr = j). We expect 
ci(oo) = Ci and a{oo) ~ 8. FIG. 3 provides an example 
of such a fit with ^ = 40, showing excellent agreement 
between numerical results and the Bcsscl function 
form of the correction to saturation. FIG. 4 shows 
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FIG. 3: (color online) Correction to the saturation of the en- 
tanglement entropy of the gapped Ising chain for ^ = 40. The 
solid circles represent the numerical values after subtraction 
of the exact saturation terms. The dashed red line repre- 
sents the function |i5'o(2L^~^). The agreement is striking 
and clearly improves for increasing values of L^~^. 

the data for ci(^). The points agree very well with 
the exact function of ^ predicted in [16, 19], lending 
support to our method. Finally, FIG. 5 illustrates the 
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FIG. 4: (color online) The circles are numerical values of ci (^) 
obtained by fitting (6). The dashed line represents exactly 
twice the function in equation (13) of [19], as expected since 
Peschel's work considers a partition into two semi-infinite re- 
gions. The agreement is again extremely good. 

main result for the Ising chain. We show the function 
a(^) =: ax{^) obtained from a fit of (6) in regions 
L > x£, ioi X = 1.5,2,2.5,3,3.5,4 and 4.5. In all cases, 
CfxiO appears constant as a function of ^, and as x 
increases this constant approaches a{oo) — 8. 
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FIG. 5: (color online) Numerical values of Qi(^). Linear fit- 
tings of the seven sets of data provide the asymptotic values 
ai.s{oo) = 7.85581, 02(00) = 7.9186, 02.5(00) = 7.95252, 
03(00) = 7.97121, 03.5(00) = 7.9818, 04(00) = 7.98783 and 
04.5(00) — 7.99146. Extrapolating in x gives the prediction 
lim,^ooO:,(oo) = 7.99281. 



XXZ spin 1/2 model in the gapped regime. Let 

us now consider the Hamiltonian (with J > 0) 

N 

i/ = J ^ (af a?+i + af a^i + Aa^a^i) , (7) 



with periodic boundary conditions. This is the anti- 
ferromagnetic spin-i XXZ chain (anisotropic Heisenberg 
model). This model displays a rich variety of features, 
depending on the value of the anisotropy parameter A. 
In the thermodynamic limit N ^ 00 the model has a 
critical line on A G [—1,1]. For A G (—1,1], includ- 
ing the Heisenberg model at A = 1, the thermodynamic 
limit is well described by a CFT (a free massless bo- 
son) with c = 1. The other critical point A = —1 is 
not described by CFT but rather corresponds to a crit- 
ical theory with an infinitely degenerate ground state; 
the entanglement of such states has been studied in [20] . 
For A > 1 the model is gapped, hence the scaling limit 
should be described by a massive QFT. The QFT is the 
S'C/(2)-Thirring model (the sine-Gordon model at a spe- 
cial value of its coupling constant - different values of the 
coupling are recovered by approaching the critical line in 
other ways, see e.g. [21, 22]). This model has a spectrum 
of two asymptotic particles of equal mass. 

In the XXZ model the exact correlation length is 
given by ^{A)'^ = 2 + E^i ^"^ tanh(n7), where 
7 = cosh"^(A) [18]. 

Numerical simulations on the XXZ chain were obtained 
by enrploying the density matrix renormalization group 
(DMRG) approach [24] . In contrast to the Ising example, 
a new set of technical challenges arises. Whereas in the 
Ising model we could use formulae for the entropy where 
the length of the chain was infinity from the outset, in 
the XXZ case we have to deal with finite chains. This 
means that if we want to get meaningful results in the 
scaling limit we need to ensure that N :^ L and, at the 
same time, if we want to test the behaviour (3) we have 
to consider L ^ ^ whilst ^ ^ e is sufficiently large so as 
to be close to the critical point A = 1. 

We have found that for a given A an optimal choice is 
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FIG. 6: (color online) The function Ci(^) near the critical 
point. The points represent the numerical results; the dotted 
line, the exact function as predicted in [22, 23]. Denoting by 
f/(5) the function (14) of [22] minus the leading logarithmic 
term |log^, the dotted line above is the function ci(^) = 
2C/(^) + log 2. The factor 2 is due to the presence of two 
boundary points. The term log 2 is related to the degeneracy: 
our zero-momentum ground state is composed in equal parts 
of two Neel-like orthogonal states of the same entropy. This 
function has the property ci(co) = 2 log 2/3. 



For 1/ ^ ^ we observe small oscillations in the en- 
tropy between blocks with odd and even numbers of 
sites. These oscillations are slightly visible in the co- 
cfhcient ci(^), FIG. 6, and the Bessel function, FIG. 7, 
and clearly visible in the function a(^), FIG. 8. Oscilla- 
tions in the entropy have been seen for finite chains (see 
e.g. [26]) where they arise as a consequence of finite size 
effects, albeit not for the von Neumann entropy. Clearly, 
in our case these oscillations must have a different origin. 
We believe that they are a consequence of the degener- 
acy of the ground state: for A > 1 the ground state is 
two-fold degenerate [18, 25] in the thermodynamic limit 
(diagonalizing the momentum, two eigenvalues and tt 
occur). Since we are always dealing with finite chains 
we should never see this degeneracy in our numerical 
simulations (the state with momentum eigenvalue is 
the "true" ground state). However for long chains, the 
gap between the two states narrows and the precision of 
our algorithm fails to distinguish the energies of those 
states. The oscillations should then be due to an al- 
ternate targeting of two linear combinations of the true 
(0-momentum) ground state and the 7r-momentum state 
(with small amplitude), depending on the parity of L. 



achieved by setting N = 5^. Longer chains are too hard 
to simulate whilst for shorter chains boundary effects 
make it impossible to obtain meaningful results. We con- 
sidered the cases ^ = 12, 14, 16, 18, 20 and 22, and used 
again the form (6). We fitted over the range £_ < L < 2£_. 

A further challenge is posed by the fact that we have 
to consider periodic boundary conditions to recover the 
behaviour (6). This makes the convergence of our DMRG 
algorithm much slower, forcing us to consider up to 1000 
states to observe good convergence. 

In order to guarantee the correct identification of the 
ground state we have employed two control parameters: 
the local magnetization (af) and the truncation error. 
We used the condition (erf) = as guiding principle for 
selecting the ground state. Imposing this condition is 
quite challenging for periodic chains. In our work we 
have been able to guarantee (erf) < 10^-^. 
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FIG. 7: (color online) Correction to the saturation of the en- 
tanglement entropy of the gapped XXZ chain for ^ = 20. The 
solid circles represent the numerical values after subtraction 
of the exact saturation terms. The dashed red line represents 
the function jKo{2L^~^). The agreement is good and clearly 
improves for increasing values of L^~^. 
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FIG. 8: (color online) Values of a(^) approaching the critical 
point. The linear expression q(^) = -4.3261^"^ + 4.0039 
reasonably fits the points and its extrapolation gives a(cxj) = 
4.0039, which is quite close to the QFT prediction a^^'^ = 4. 



FIG. 8 shows that the coefficient a(^) converges with 
great precision to its theoretical value 4. Contrary to 
the Ising case, a{£,) has a dependence on the correlation 
length. FIG. 7 shows a good fit to a Bessel function. 
Interestingly, this fit continues to be good away from A = 
1 as long as the coefficient a(^) is changed as dictated by 
FIG. 8. This suggests that the function (6) may provide 
to a good accuracy a universal description of the entropy 
of gapped spin chains. 

Finally, at the critical point, we have found that the 
entropy is very well fitted by the function (2) with c ~ 
1.00024 and C2 = 0.733758. We can therefore make a 
prediction for the universal QFT constant U = ci{oo) — 
C2 = -0.27166 in the S'C/(2)-Thirring model. 

Conclusions. By studying gapped spin chains we 
have provided strong numerical evidence for the be- 
haviour (3) suggested in [1, 2] for two-dimensional QFT. 
Besides confirming QFT predictions, our results show 
that computing the entanglement entropy may be a good 



numerical way to determine the number of lightest parti- 
cles in the universal spectrum, especially in combination 
with the knowledge that a(oo) in (6) must be a fraction 
of 8. Our results also suggest that some of the QFT pre- 
dictions may still hold, with small changes, beyond the 
scaling regime. It would be very interesting to derive this, 
or the appropriate modification of (6), for spin chains 
from first principles. We have made a prediction for the 
universal constant U in the S'[/(2)-Thirring model, still 
to be confirmed by means of QFT techniques. Finally, it 
would be very interesting to carry out a similar analysis 
for non-integrable spin chains and thus verify the main 
result of [2], namely, the fact that the behaviour (3) also 
holds for 1-1-1 dimensional non-intcgrablc QFTs. 
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